Objetivo: aprender cómo generar tus propias funciones y realizar cálculos paralelos.
R permite al usuario crear funciones. El modo de estas expresiones es function, se almacenan en una forma especial y se pueden utilizar en expresiones posteriores.
Una función se define de la siguiente manera:
nombre_funcion nombre de la función definida por el usuario.
argn son los diferentes argumentos de la función.
expresion es el proceso que seguirá la función.
La función se aplica de la siguiente manera:
Como ejemplo, generemos una función que cuente los caracteres en una cadena de texto. Primero comenzaremos generando el esqueleto de la función y agregando un parámetro llamado word:
Luego, aplicamos la función nchar.
Finalmente, podemos incluir un mensaje con el resultado. Observa que utilizaremos la función return para indicar explícitamente la salida de nuestra función.
Ahora apliquemos la función:
count_char(word = "Oscar Baez")
## [1] "La palabra Oscar Baez tiene 10 caracteres!"
count_char(word = "Análisis espacial")
## [1] "La palabra Análisis espacial tiene 17 caracteres!"
count_char(word = 10:14)
## [1] "La palabra 10 tiene 2 caracteres!" "La palabra 11 tiene 2 caracteres!"
## [3] "La palabra 12 tiene 2 caracteres!" "La palabra 13 tiene 2 caracteres!"
## [5] "La palabra 14 tiene 2 caracteres!"En el último ejemplo, proporcionamos un vector numérico al parámetro word. Podríamos controlar la entrada del parámetro. Por ejemplo:
count_char <- function(word){
if(!is.character(word))
stop("¡La entrada no es un caracter!")
n <- nchar(word)
n <- paste("La palabra", word, "tiene", n, "caracteres!")
return(n)
}La función stop interrumpirá la ejecución de la función si no se cumple la condición.
¿Cómo crearías una función que te devuelva el valor máximo en un vector dado?
Si lo aplicamos:
Ahora que hemos creado nuestra propia función, ¿cómo podemos aplicarla a un conjunto de datos raster para obtener el valor máximo por celda de grilla durante un período definido? Para este propósito, podemos utilizar la función app del paquete terra.
Para este propósito, utilizaremos datos mensuales CHIRPSv2 para el año 2000.
print(chirps)
## class : SpatRaster
## dimensions : 2000, 7200, 12 (nrow, ncol, nlyr)
## resolution : 0.05, 0.05 (x, y)
## extent : -180, 180, -50, 50 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## sources : chirps_monthly2000-01.tif
## chirps_monthly2000-02.tif
## chirps_monthly2000-03.tif
## ... and 9 more source(s)
## names : chirp~00-01, chirp~00-02, chirp~00-03, chirp~00-04, chirp~00-05, chirp~00-06, ...
## min values : 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, ...
## max values : 1251.456, 1467.558, 1222.492, 1224.301, 1322.237, 1584.747, ...Ahora podemos aplicar la función app para obtener el valor máximo en cada celda.
Calcularemos la precipitación media anual utilizando 12 archivos mensuales de CHIRPSv2 sobre los países de la cuenca del Nilo en este ejemplo:
print(countries)
## class : SpatVector
## geometry : polygons
## dimensions : 15, 9 (geometries, attributes)
## extent : 12.19545, 47.98618, -13.45568, 31.66734 (xmin, xmax, ymin, ymax)
## source : Nile_Basin_Countries_GAUL2014_2.shp
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## names : STATUS DISP_AREA ADM0_CODE ADM0_NAME STR0_YEAR EXP0_YEAR
## type : <chr> <chr> <int> <chr> <int> <int>
## values : Member State NO 205 Rwanda 1000 3000
## Member State NO 133 Kenya 1000 3000
## Member State NO 74 South Sudan 2011 3000
## Shape_Leng Shape_Area Name_label
## <num> <num> <chr>
## 8.127 2.064 Rwanda
## 48.55 47.35 Kenya
## 46.91 51.6 South SudanPor lo tanto, queremos crear una función que:
Una vez que la función esté lista, se puede aplicar a las diferentes entidades del shapefile.
Primero, podemos preparar los datos que se utilizarán como entrada en la función. Guardaremos las rutas de los rasters en un objeto llamado raster_paths y el shapefile en un objeto llamado shape.
print(shape)
## class : SpatVector
## geometry : polygons
## dimensions : 1, 9 (geometries, attributes)
## extent : 28.86187, 30.89965, -2.839881, -1.047913 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## names : STATUS DISP_AREA ADM0_CODE ADM0_NAME STR0_YEAR EXP0_YEAR
## type : <chr> <chr> <int> <chr> <int> <int>
## values : Member State NO 205 Rwanda 1000 3000
## Shape_Leng Shape_Area Name_label
## <num> <num> <chr>
## 8.127 2.064 RwandaLuego, podemos crear la función. La llamaremos get_map para la precipitación media anual.
Esta función requiere dos entradas: las 12 rutas de los rasters globales de precipitación (uno para cada mes del año) y el shapefile. Podemos proporcionar las rutas de los objetos y leerlos dentro de la función.
Al final, asignaremos estos objetos a la función, es decir, raster_paths = raster_paths y shape = shape. El primer paso sería importar los datos raster utilizando la función rast del paquete terra.
El segundo paso sería recortar y enmascarar la pila de rasters a la extensión del shapefile.
Finalmente, podemos sumar las 12 capas de precipitación y devolver el resultado.
Una vez que tenemos la función lista, podemos aplicarla:
El procesamiento de grandes cantidades de datos espaciales puede llevar mucho tiempo. Sin embargo, la mayoría del código de R se ejecuta en un solo procesador. La mayor parte del tiempo esto no es un problema, pero a veces estod procesos pueden:
Requerir bastante tiempo
Consumir memoria
Tomar bastante tiempo al leer o escribir archivos
Tomar bastante tiempo al transferir datos
Las computadoras tradicionales tienen un solo CPU, que a su vez puede contener múltiples núcleos. Estos procesadores y núcleos pueden realizar cálculos (ten en cuenta que las computadoras modernas pueden tener múltiples procesadores).
La siguiente función genera una lista de números primos hasta un número determinado. Esta función fue tomada de stackoverflow.
prime_numbers <- function(n){
n <- as.integer(n)
if(n > 1e6) stop("n too large")
primes <- rep(TRUE, n)
primes[1] <- FALSE
last.prime <- 2L
for(i in last.prime:floor(sqrt(n)))
{
primes[seq.int(2L*last.prime, n, last.prime)] <- FALSE
last.prime <- last.prime + min(which(primes[(last.prime+1):n]))
}
which(primes)
}
prime_numbers(100)
## [1] 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97Utilizaremos la función Sys.time para evaluar el tiempo necesario que lleva esta función para obtener números primos para una secuencia que comienza en 10 y termina en 15,000.
Para paralelizar el código, utilizaremos el paquete doParallel. Este paquete proporciona funciones para la ejecución paralela de código R en máquinas con múltiples núcleos o procesadores o múltiples computadoras.
Tenemos que indicar el número de núcleos que se utilizarán. Para ello, podemos usar la función detectCores.
No se recomienda utilizar todos los núcleos para que tu computadora no se bloquee. Por lo tanto:
Después, tenemos que registrar el backend paralelo. Para ello, utilizaremos la función registerDoParallel.
Finalmente, el proceso iterativo del for for será diferente. En lugar de for, utilizaremos foreach; en lugar de in, utilizaremos =. A continuación, podemos indicar explícitamente qué paquetes se utilizarán con el parámetro .packages (repasaremos esto en el ejemplo). Finalmente, el operador %dopar% debe colocarse antes de abrir el foreach.
La función stopImplicitCluster se puede utilizar en documentos y otros lugares donde es importante cerrar explícitamente el clúster creado implícitamente.
Retomando el ejemplo de extracción de precipitación media anual por país, obtengamos la precipitación media anual para todos los países de la cuenca del Nilo.
Al observar los datos, podemos ver que hay áreas en disputa en el shapefile. Como solo queremos los países, podemos excluirlas.
countries <- countries[which(countries$DISP_AREA == "NO"),]
data.frame(countries)
## STATUS DISP_AREA ADM0_CODE ADM0_NAME STR0_YEAR
## 1 Member State NO 205 Rwanda 1000
## 2 Member State NO 133 Kenya 1000
## 3 Member State NO 74 South Sudan 2011
## 4 Member State NO 257 United Republic of Tanzania 1000
## 5 Member State NO 253 Uganda 1000
## 6 Member State NO 79 Ethiopia 1000
## 7 Member State NO 77 Eritrea 1000
## 8 Member State NO 43 Burundi 1000
## 9 Member State NO 68 Democratic Republic of the Congo 1000
## 10 Member State NO 6 Sudan 2011
## 11 Member State NO 40765 Egypt 1000
## EXP0_YEAR Shape_Leng Shape_Area Name_label
## 1 3000 8.127003 2.063852 Rwanda
## 2 3000 48.549430 47.345770 Kenya
## 3 3000 46.905431 51.599166 South Sudan
## 4 3000 82.950601 76.860266 United Republic of Tanzania
## 5 3000 23.244416 19.629674 Uganda
## 6 3000 50.380131 92.869258 Ethiopia
## 7 3000 50.005783 10.167958 Eritrea
## 8 3000 9.118459 2.185652 Burundi
## 9 3000 98.951177 189.998107 Democratic Republic of the Congo
## 10 3000 81.910242 155.888802 Sudan
## 11 3000 61.251157 89.079113 EgyptAhora, podemos aplicar la función get_map para cada país. Dado que este proceso seguramente llevará mucho tiempo, podemos paralelizarlo. Primero, aquí está la función get_map.
get_map <- function(raster_paths, shape){
# leyendo los datos
p_rast <- terra::rast(raster_paths)
# recortando y enmascarando el conjunto de datos ráster
p_rast <- terra::crop(p_rast, shape, snap = "out")
p_rast <- terra::mask(p_rast, shape)
# sumando las capas de precipitación
p_rast <- sum(p_rast)
return(p_rast)
}Primero, almacenamos las rutas de los archivos ráster en un objeto como hicimos antes. Luego, podemos establecer la ruta de salida donde queremos almacenar los archivos resultantes.
Construimos el bucle foreach de la siguiente manera:
cores <- detectCores() - 1
registerDoParallel(cores = cores)
foreach(i = 1:nrow(countries), .packages = "terra") %dopar% {
r <- get_map(raster_paths = raster_paths, shape = countries[i,])
n <- countries[i,]$Name_label
n <- paste0(output, "/", n, ".tif")
writeRaster(r, n, overwrite = TRUE)
}
stopImplicitCluster()Después del proceso, deberíamos tener los rásteres en la carpeta de salida.
